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Abstract 

A model consisting of a mixed Kuramoto - Sivashinsky - Korteweg - de Vries 
equation, linearly coupled to an extra linear dissipative equation, is proposed. 
The model applies to the description of surface waves on multilayered liquid 
films. The extra equation makes its possible to stabilize the zero solution in the 
model, thus opening way to the existence of stable solitary pulses (SPs). By 
means of the perturbation theory, treating dissipation and instability-generating 
gain in the model (but not the linear coupling between the two waves) as small 
perturbations, and making use of the balance equation for the net momentum, 
we demonstrate that the perturbations may select two steady-state solitons from 
their continuous family existing in the absence of the dissipation and gain. In 
this case, the selected pulse with the larger value of the amplitude is expected to 
be stable, provided that the zero solution is stable. The prediction is completely 
confirmed by direct simulations. If the integration domain is not very large, some 
pulses are stable even when the zero background is unstable. An explanation 
to the latter finding is proposed. Furthermore, stable bound states of two and 
three pulses are found numerically. 

PACS: 05.45.Yv, 46.15.Ff 
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1 Introduction 



The Kuramoto - Sivashinsky (KS) equation is a well-known model of one- 
dimensional turbulence, which was derived in various physical contexts, includ- 
ing chemical-reaction waves [Q, propagation of combustion fronts in gases |2|, 
surface waves in a film of a viscous liquid flowing along an inclined plane |3|, 
patterns in thermal convection Q, rapid solidification ||^, and others. It has 
the form 

Ut + UUx = -aUxx - lUxxxx, (1) 

where a > and 7 > are coefficients accounting for the long-wave instability 
(gain) and short-wave dissipation, respectively. 

A generalized form of the KS equation contains a linear dispersive term 
borrowed from the Korteweg - de Vries (KdV) equation, 

~t- UUx ~t~ ^xxx — ^"^xx ^^xxxx • (2) 

As well as the KS equation proper, the generalized equation (||) applies to the 
description of surface waves on flowing liquid films |^, and it also serves as 
a general model which allows one to study various nonlinear dissipative waves 
particular, a subject of considerable interest was the study of solitary- 
pulse (SP) solutions to both the KS equation Q and Eq. (||) 10 . By means 



of numerical methods, it is possible to find a vast family of SP solutions to Eqs. 

and in the form u{x^ t) = u(x — st) with a constant velocity s. However, 
it is obvious that all these solutions are unstable in an infinitely long system, as 
the zero solution, into which SP goes over at \x — st\ — > 00, is unstable in both 
equations. 

It is an issue of a principal interest to find a physically relevant model which 
combines the dissipative and dispersive features, and simultaneously supports 
stable SPs. It appears that the simplest possibility to construct such a model is 
to couple Eq. (H) to an extra linear stabilizing equation, arriving at a system 

Ut + UUx + Uxxx = -aUxx - "fUxxxx + f-lVx, (3) 
Vt + CVx = TVxx + ^2Ux, (4) 

where the dissipative parameter (effective diffusion coefficient) F > accounts 
for the stabilization (see below) , and c is a group- velocity mismatch between the 
two wave modes. The coupling parameters ei and £2 must have the same sign 
(otherwise the coupling gives rise to an instability), while their magnitudes may 
be different. However, it is always possible to make them equal, ei = £2 = e, 
by means of an obvious rescaling of u and v. Then, using the remaining scaling 
invariance of the equations, it is possible to set e = 1. Thus, we will be dealing 
with a system containing four irreducible parameters, 

Ut + UUx + Uxxx -Vx = -OiUxx - lUxxxx, (5) 
Vt+CVx-Ux = TVxx- (6) 
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Note that the system conserves two "masses" , 

/+CX; p + co 

u{x)dx, N= v{x)dx. (7) 

-oo J —oo 

The system of equations (H) and can find their natural physical realiza- 
tion as a model describing coupled surface and interface waves in a two-layered 
flowing liquid film, cf. the similar interpretation of the single equation (|l]) or Eq. 

mentioned above. In particular, the linear coupling via the first derivatives 
is the same as in known models of coupled internal waves propagating in multi- 
layered fluids Then, the linear dissipative equation (^) implies that the 
substrate layer is essentially more viscous than the upper one. In fact, the addi- 
tional equation may also be nonlinear, but it can be checked that inclusion 
of the nonlinear term vvx into this equation docs not produce any conspicuous 
difference, therefore we focus on the simplest model (||), (||), which provides for 
the stabilization of SPs. 

The system of equations (||) and (||) is qualitatively similar to a system of 
linearly coupled Ginzburg-Landau (GL) equations describing propagation of lo- 
calized pulses in an fiber-optic core equipped with distributed gain, which is 
linearly coupled to an extra lossy core that provides for the stability of the 
pulses H, 0, |l| (such double -core systems have recently become available to 
experimental studies, and they have very promising features for applications to 
optical communications, see a short overview in Ref. p^). The most funda- 
mental version of this GL system is that in which the extra stabilizing equation 
is also linear, cf. Eq. (||); in that case, SP solutions can be found in an exact 
analytical form, and they are stable in a certain parametric region p^ . 

In this work, we will find stable SPs in the system (|), (|), which appears 
to be the first example of stable pulses in a model of the KS type. In section 2, 
we analyze the stability of the zero solution, which, as it was mentioned above, 
is a necessary condition for the stability of SPs in the infinitely long system. In 
section 3, an analytical perturbation theory for the pulses is developed, which 
is based on treating the gain and dissipation constants a, 7, and F in Eqs. 
and (|^) as small parameters (while the group-velocity mismatch c need not be 
small). In the zero-order approximation, a = 7 = F = 0, Eqs. (|^) and (^ have 
a one-parameter family of exact soliton solutions. Using the known approach 
based on the balance equation for the momentum |7|, we demonstrate that the 
combination of the perturbation terms in Eqs. (|) and (^ may select one or 
two stationary pulses out of the continuous family existing in the zero-order 
approximation. As is known [Q, the existence of two different SP solutions 
is a necessary condition for the stability of one of them, the second pulse (the 
one with smaller amplitude) being unstable, as it plays the role of a separatrix 
between attraction domain of the zero solution and stable pulse. 

In section 4, we present results of direct numerical simulations of the full 
system (||), (|6|), which demonstrate that stable SPs exist indeed. In fact, simu- 
lations sometimes produce stable pulses even in the case when the zero solution 
is not stable. This stability extension may be explained by the finite size of 



3 



the simulation domain. Moreover, stable bound states (BSs) of two and three 
pulses are found numerically too. 



2 The stability of the zero solution 

As it was explained above, it is necessary to investigate the stability of the 
trivial solution, u = v = 0, before the consideration of pulses. To this end, 
we substitute into the linearized equations (||) and (^) a perturbation in the 
form u ~ exp {ikx + Xt), v ~ exp {ikx + \t), where k is an arbitrary real wave 
number of the perturbation, and A is the corresponding instability growth rate, 
which leads to a dispersion equation, 

(A - ik^ ~ ak^ + ^k^) (A + ic± + rfc2) + fc2 = 0. (8) 

The stability condition states that both solutions of the quadratic equation (||) 
must satisfy the inequality Re [A(fc)] < at all the real values of k. 

For fc — > 0, solutions to Eq. (ph can be found in the form of an expansion 

\{k) ^ iXik + X2k^ + ... , (9) 

where Ai = (— c ± Vc^ + 4) /2, and 



- ' ^ 

the sign ± being the same in Ai and A2. As Ai is always real, at the first 
order the expansion (|^) implies neutral stability. The expression (10) yields a 
necessary condition for the stability of the zero solution at the second order of 
the expansion. Re A2 < 0, which can be cast into a form 

r-a>\/af|c|. (11) 

In the particular case c = 0, the condition ( pT| ) amounts to F > a, which has 
a simple meaning: the stabilizing diffusion coefficient in Eq. (H) must be larger 
than the instability-driving "anti-diffusion" (gain) coefficient in Eq. (||) . A very 
similar necessary stability condition is known in the above-mentioned system 
of coupled GL equation describing a dual-core optical fiber with one active and 
one passive cores |l4| . 

Comprehensive analysis of the zero-solution stability was performed by means 
of numerical solution of the dispersion equation (||) . It was found that the full 



stability condition does not amount to the inequality (11) (i.e., it may happen 
that Re [A(fc)] is negative at small fc, but it takes positive values in some inter- 
val of finite values of k). Numerically found stability borders in the plane of 
the parameters (a,F) for a fixed value 7 = 0.05 of the short-wave stabilization 
parameter in Eq. (^ and two different values of the group-velocity mismatch, 
c = and c = — 1, are shown by dashed curves in stability diagrams for the 
pulses displayed in Figs. 1 and 2. In both cases, the zero solution is stable to 
the left of the stability border. Note that, for small a and F, the zero-solution 
stability region is indeed determined by Eq. (^ , but at larger values of a and 
F there appear additional stability restrictions. 
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3 The perturbation theory for sohtary pulses 



At the zeroth order, setting 7 = F = a = in Eqs. (||) and but keeping an 
arbitrary value of c, we arrive at a conservative system consisting of the KdV 
equation coupled to an extra linear one, 

Ut + UUx + Uxxx = Vx,Vt + CVx ^ u^. (12) 

Equations (|2|) have a family of exact two-component soliton solutions, 

u{x, t) — 1277^sech^ (77 {x — st)) , v{x, t) — [c — s) ^ ■ u{x, t), (13) 

where 77 is an arbitrary parameter which determines the soliton's amplitude and 
width, and the velocity s takes two different values for given 77, 



(c + V) =F y^(c- 4772)2 + 4 



(14) 



It will be more convenient to use, as a parameter of the soliton family, not the 
amplitude 7/, but rather the relative velocity, 

S = c-s, (15) 

in terms of which the amplitude is given by an expression obtained from Eq. 

477^ =c-(5 + l/(5. (16) 

The range of meaningful values of S is restricted by the condition 77^ > 0. 

We have checked by direct simulations of Eqs. (|lj) that all the soliton solu- 
tions (|l^) are stable within the framework of the unperturbed equations (|l^ . 
On the other hand, simulations also clearly demonstrate that collisions between 
solitons having different velocities are inelastic (although not strongly inelastic, 
see a typical example in Fig. 3), hence the conservative system (p^, unlike the 
KdV equation proper, is not an exactly integrable one. The nonintegrability of 



the system (12) has also been confirmed by analysis of its symmetries performed 
by G. Burde (unpublished). 

The next step is to restore the small dissipative and gain perturbations, 
getting back from Eqs. ( p^ ) to Eqs. (^) and (|6|). To this end, we notice that 
the unperturbed equations (|l^) conserve not only the masses (^ but also the 



net momentum, 

,2 I „,2 



P = i / (u^+ v^) dx. (17) 
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Following the work , in the first approximation of the perturbation theory the 
evolution of the soliton may be described by means of the balance equation for 
the momentum. Indeed, a consequence of Eqs. (^ and (||) is the following exact 
evolution equation for the net momentum in the presence of the perturbations: 



dP_ 
'dt 



-\-oc 

{aul - jul^ - Tvl) dx. (18) 
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Steady-state SPs are selected by the condition dP/dt = 0. The right-hand 
side of Eq. (|l|) can be exphcitly calculated in the approximation in which u 
and V are substituted by the expressions (|l^). After a straightforward algebra, 
the equation dP/dt = can be cast into the form of a cubic equation for the 
relative velocity 6 of the unperturbed soliton (|l^) , 

5(5^ -I- (75 - 5c) ^2 _ 5J _ 7f = 0, 5 = a/7, f = T/j. (19) 

Roots of Eq. (H) select SPs that may exist as steady states within the 
framework of the perturbation theory. Note that, besides the obvious condition 
that physical roots for 6 must be real (they may be both positive and negative), 
they must also satisfy a condition that, after the substitution into Eq. (p^, 
they must produce 77^ > 0. Generally speaking, there may exist up to three 
physical roots of Eq. (p^); however, in a vast parametric area considered, we 
have never encountered a case when Eq. ( pj| ) would indeed have three physical 
roots, while the existence of two physical solutions is quite possible, see below. 
As it was mentioned in the Introduction, one may expect that SP may be stable 
if precisely two different pulses exist, then the one with the larger amplitude ry^ 
has a chance to be stable, while the pulse with the smaller amplitude is always 



unstable 14 . Indeed, if there is a stable SP, we are dealing with a bistable 
system, as the parameters are chosen so that the zero solution is also stable, 
see the previous section. In a bistable system, there should exist a separatrix, 
i.e., a border between attraction domains of two stable solutions, the separatrix 
itself being an unstable stationary solution. In the situation with two different 
stationary SP solutions predicted by the perturbation theory, the one with the 
smaller amplitude and larger width is a natural candidate to the role of the 
unstable separatrix solution, while its counterpart with the larger amplitude 
and smaller width may be stable (generally speaking, it may be stable in a part 
of the parametric region where this situation takes place ||l^, |l^ ) . 

Equation (|l^) may be simplified in the case c = 0, if we additionally assume 
that both renormalized parameters a and F are large (in fact, we are interested 
in the case when T ^ 3> 1). Then, the term —5S may be neglected in Eq. 
^§1), so that it takes the form 

56^ + 75(5^ - 7r = 0. (20) 

In the case c = 0, the condition rj^ > following from Eq. (^6|) also takes a 
simple form: a physical root is that which belongs to either of the two intervals, 

(5 < -1; < (5 < 1. (21) 

With regard to the assumption that F and a are large, it is easy to see that 
the simplified equation ( pO| ) always has a real root in the region S > 1, which is 
unphysical according to Eq. (pT]). Two physical roots (5 < — 1 exist under the 
condition 
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Note that this condition does not contradict the necessary condition (^l|) of the 
stabiHty of the zero solution, which in the present case (c — 0) takes the form 
r > 3. Indeed, the latter inequality is compatible with (|2|), provided that 
5 > -\/3 (15/14) sa 1. 856, which is correct, as we are here dealing with the case 
when a is large. However, the inequality ( p^ ) is not necessarily compatible with 
the full stability condition for the zero solution, see Fig. 1. 

In the general case, it is easy to solve Eq. numerically. Then, selecting 
a parametric region in which there are exactly two physical solutions (which, as 
it was explained above, is a necessary condition for the existence of a stable SP), 
one may identify a narrower region in which this condition holds and, simul- 
taneously, the zero solution is stable. Stable pulses may exist only inside that 
region where both necessary stability conditions overlap, and direct simulations 
show that all the pulses belonging to the region are stable indeed, at least in 
case displayed in Fig. 1, see details below. 

The so defined regions in the parametric plane (a,T), in which stable SPs 
are expected, are displayed, for 7 = 0.05, in Figs. 1 and 2 for c = and c = — 1, 
respectively. The condition of the existence of exactly two different physical 
solutions for the pulses holds to the right of the continuous curve in these figures 
[note that the part of the curve corresponding to sufficiently large values of F 
is well approximated by the analytical expression (^) obtained above] . At the 
points belonging to the continuous curve, the two physical solutions merge and 
disappear via a typical tangent (saddle-node) bifurcation. 

The same analysis performed for values of the short-wavelength-dissipation 
parameter 7 different from the value 0.05, for which the results are presented 
in Figs. 1 and 2, shows that variation of 7 produces little change in terms of 
the expected SP stability region (generally, the size of the region increases with 
7) . As for the effect of the group- velocity mismatch c, we have found that the 
stability region quickly shrinks with the increase of c when c is positive, and 
there is no stability region at c > Ccr, where Ccr is slightly larger than 0.3. In 
that case, the areas in which, respectively, the zero solution is stable, and there 
are two different stationary SPs, do not overlap. To illustrate this point, a very 
narrow stability region in the (a,r) plane for c — 0.3 is shown in Fig. 4. 

4 Numerical simulations of the solitary pulses 

As it was said above, it is necessary to directly check whether stable SPs indeed 
exist in the region where the stability is expected. To this end, Eqs. (^), (||) 
with periodic boundary conditions (b.c.) were integrated by means of an implicit 
Fourier spectral method ||l^, the time step being, typically, 0.01 and 0.02. (a 
description of the method is given in the appendix) . The initial conditions were 
taken as suggested by the perturbation theory, i.e., in the form (p^), but with 
arbitrary values of the amplitude, in order to check whether strongly perturbed 
pulses relax to stable ones, i.e., whether the stable pulses are attractors. 

Results are displayed in Fig. 1 by means of the symbols x and o, standing for 
unstable and stable solitons, respectively (it may happen that, near the border 
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with the unstable SPs, some pulses which appear to be stable are subject to 
a very weak instability which does not manifest itself within the integration 
time limits). As is seen, the pulses are indeed stable everywhere inside the 
expected stability region. Moreover, all the stable pulses were found to be 
strong attractors. For instance, in the case a = 0.1, F = 0.15, the initial pulses 
definitely relaxed to a single stationary SP if their initial amplitude Aq exceeded 
1.7. In particular, starting with = 3 and Aq = 12 at t — 0, the pulse attained 
the value of the amplitude, respectively, A = 6.91 and A ~ 7.18 at t — 400. The 
analytical prediction for the amplitude of the steady-state pulse (13) (the one 
with the larger value of the amplitude) is, at the same values of the parameters, 
^anai = 121]^ « 6.45, SO that a discrepancy with the numerical results is less than 
10%. On the other hand, if the initial amplitude was too small, e.g., Aq — 0.75, 
the pulse decayed into zero, which is natural too, as the zero solution has its own 
attraction basin. Note that for the second (smaller) steady-state pulse, which 
is expected to play the role of the separatrix between the attraction basins of 
the stable pulse and zero solution, the perturbation theory predicts, in the same 
case, the amplitude Aanai ~ 2.15, so it seems quite natural that the initial pulses 
with Aq ~ 3 and = 0.75 relax, respectively, to the stable pulse and to zero. 

Figure 1 shows that the numerically found upper border of the stable- 
pulse region is quite close to the border of the existence region for the steady- 
state pulses, as predicted by the perturbation theory. Unlike this, the nu- 
merically identified stability region extends far below the analytically found 
border of the zero-solution instability. For instance, it was found that, at 
a = 0.15 and F — 0.2, when the zero solution is unstable against perturbations 
with finite wavenumbers fc, the fastest growing perturbation corresponding to 
k = fcmax ~ 1.3, a fairly stable pulse with the amplitude A = 11.75 was found 
in the simulations, the pulse's amplitude predicted by the perturbation theory 
being j4anai ~ 11.61 in this case. Moreover, we ran simulations in which the 
most dangerous perturbation with the above-mentioned wavenumber A: = 1.3 
and a rather large amplitude, Aport = 1, was deliberately added to the pulse 
in the initial configuration. Instead of growing and destroying the pulse, the 
perturbation was suppressed by the pulse, which remained stable indefinitely 
long (Fig. 5). 

A cause of this extended stability can be understood. A similar feature was 
reported, for a generalized asymmetric (with respect to the reflection x —x) 
cubic-quintic GL equation with periodic b.c, which has moving-pulse solutions, 
in Ref. ||l^. An explanation was that the moving pulse, traveling across the 
integration domain, periodically passes each point and suppresses the perturba- 
tion at a rate which exceeds the perturbation growth rate (see also Ref. JTof , 
where stable pulses were observed in the KS - KdV equation; in that work, an 
explanation was that the moving pulse was able to escape growing perturbation 
wave packets). It seems very plausible that a similar "sweeping" mechanism 
explains the anomalous pulse stability in the present model. Indeed, when we 
repeated the simulations for the same values of the parameters but in a spatial 
domain four times as large (i.e., the corresponding sweeping rate is four times 
as small), the pulse demonstrated the expected instability, even without any 
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specially added perturbation seed, see Fig. 6. 

Similar to a result reported for a system of linearly coupled Ginzburg-Landau 
(GL) equations in Ref. bound states (BSs) with two or more peaks can be 
found in the present model, in addition to single SPs. To this end, we took an 
initial configuration constructed as a set of two or more identical exact solutions 
of the unperturbed system ( p^ with a certain separation between them. The 
simulations have shown that BSs featuring two or more peaks of equal ampli- 
tudes indeed develop and propagate stably. These results for the two-peak and 
three-peak BSs are illustrated by Figs. 7 and 8, respectively. We have also 
checked that, even if the amplitudes of the initial pulses and separations be- 
tween them are changed, the same BS consisting of equidistant equal-amplitude 
peaks finally develops, i.e., the BSs are fairly robust dynamical objects. In this 
connection, it is relevant to mention that, in the above-mentioned coupled GL 
equations, only two-pulse BSs are completely stable, while BSs of three pulses 
are split by perturbations breaking their symmetry . 

5 Conclusion 

In this work, we have introduced a model based on the Kuramoto - Sivashinsky 
- Korteweg - de Vries equation, which is linearly coupled to an extra linear dissi- 
pative equation. The model can be applied to the description of coupled surface 
and interface waves on flowing multilayered liquid films. The additional linear 
equation makes its possible to stabilize the zero solution, which opens way to 
the existence of stable solitary pulses. Treating the dissipation and gain in the 
model (but not the linear coupling between the two wave modes) as small per- 
turbations, and making use of the balance equation for the net momentum, we 
have found that the condition of the balance between the gain and dissipation 
may select two steady-state solitons from their continuous family existing in the 
absence of the dissipation and gain (the family was found in an exact analytical 
form). When the zero solution is stable and, simultaneously, two SPs are picked 
up by the balance equation for the momentum, the pulse with the larger value 
of the amplitude is expected to be stable in the infinitely long system, while the 
other pulse must be unstable, playing the role of a separatrix between attrac- 
tion domains of the stable pulse and zero solution. These predictions have been 
completely confirmed by direct simulations. Moreover, if the integration domain 
is not very large (and periodic boundary conditions are imposed), some pulses 
are stable even when the zero background is unstable. An explanation to the 
latter feature, based on the concept of periodic suppression of the perturbations 
by the running pulse, is proposed. Furthermore, stable bound states of two and 
three identical pulses have been found numerically. An interesting issue, which 
is left for further work, is a possibility of formation of stable periodic arrays 
of the pulses. Note that periodic pulse arrays in the KS - KdV equation were 
studied in Refs. ^. 

Lastly, it is worthy to note that, since we have found stable pulses in the 
model including the dispersion term, i.e., Uxxx in Eq. (|^), a natural question is 
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if the presence of this term is a necessary condition for the existence of stable 
pulses. Our preliminary numerical study demonstrates that there is a finite 
critical minimum value of the coefficient in front of the dispersive term (if all 
the other parameters of the system are fixed) which is indeed necessary for the 
existence of stable pulses. For instance, in a typical case considered in this work, 
with c = 0, q; = 0.1,7 = 0.05, and T = 0.15, this minimum value was found to 
be roughly 0.5. 

A plausible qualitative explanation for this effect (which will be considered 
in detail elsewhere) is that, although the dispersion does not directly affect the 
saturation mechanism stabilizing the pulses, it can inhibit mode coupling as 
an effective impedance, which will lead to attaining the necessary saturation at 
a higher amplitude of the pulse. Therefore, the amplitude of pulses becomes 
larger as the dispersion increases, while the pulse widths, directly determined by 
the dissipation, remain almost constant. This trend to the formation of pulses 
with higher amplitudes in the presence of strong dispersion is a feasible cause 
for the existence of stable pulses. 
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Appendix 

Here we give a detailed description of the numerical method employed for the 
simulation of Eqs. ^ and (^. First, we introduce necessary notation for the 
discrete Fourier transform (DFT). Suppose we work on the interval / = [0,i] 
with L-periodic functions. The interval is discretized by means of a set of N 
equidistant points, Xj = jL/N, where j = 0, 1, • • • , iV — 1, with the spacing 
between them Ax — L/N. At these points, the numerical solutions for u{xj,t), 
v(xj,t) are denoted by Uj{t), Vj{t) respectively. The corresponding spectral 
variable is = ^irk/L with k e {-N/2, • • • , -1, 0, 1, • • • , N/2 - 1} (actually, 
N/2 is an integer, see below). Then DFT is given by 

, ^ , , N N 
Uk = :Fuj = 2^ u^exp{-i^kXj) , k = - — 0,1, -1. (23) 

The inverse DFT is defined as 

^ N/2-1 

Uj=T-^uu^- Uke-K^{i£,kx,), j = 0,l,---,7V-l. (24) 

k=~N/2 
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We employ the fast Fourier transform (FFT) to carry out the DFT and its 
inverse in the numerical form, therefore N must be taken equal to a powers of 
2. 

The Fourier transform converts Eqs. (^. (|^) into 

Ut + i£.k:F{^u'^)-i^lu-iS,kV = a^lu-"f£Xu, (25) 
vt + icikv - i^ku -r^fc^^: (26) 

where k = -N/2, • • • , -1, 0, 1, • • • , N/2 - 1. 



For the time integration of Eqs. (25) and (E6f) , the Crank-Nicolson scheme 



is used, which leads to a nonlinear system 

1.^ .^3lt"+l+ti" . 

77 - ^ik Kk- 



At 2 " 2 

•^(^T^)+-^(^^) =Kfe-7e^') ^ , (27) 



At 



2 

+ ic^k ;r *?fc 



-m ^ ■ (28) 



To solve this nonlinear system, the following iteration procedure is employed: 

^ ^n^ ^n+lfi ^ ^ ^29) 



i^k 7^ *Cfe- 



At ^"2 ^"2 



2 ^ \ ' 2 



£jn+l,r+l _ ^ri 



+ *c^fc *?fe- 



At ^ 2 

fn+l,r+l f fn 

= -nl- (31) 

Here r = 0, l,---,i? — 1, where R is the iteration number in each time step. 
In practice, the number of the Fourier modes was taken to be 256 or 512, the 
typical time step was 0.01 and 0.02, and the iteration was run twice in each 
time step. 
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Figure Captions 



Fig. 1. The stability region for solitary pulses in the parametric plane 
of the system (|), (|) for 7 = 0.05 and c = 0. The zero solution is stable to the 
left of the dashed curve, and Eq. (|l^) produces two physical solutions to the 
right of the continuous curve. The symbols x and o mark points at which direct 
simulations show, respectively, that the solitary pulse is unstable or stable. 

Fig. 2. The expected stability region for the solitary pulses in the parametric 
plane (a,r) for 7 = 0.05 and c = —1. The continuous and dashed curves have 
the same meaning as in Fig. 1. 

Fig. 3. A typical result of an inelastic collision between two stable solitons 
with different velocities, simulated within the framework of the zero-order con- 
servative system (||) with c = 0: (a) t = 0, (b) i = 0, (c) i = 40, (d) t = 40, 
the initial velocities of the two solitons are si = 4.236 and S2 = 2.414. 

Fig. 4. The expected nearly vanishing stability region for the solitary pulses 
in the parametric plane (a,r) for 7 = 0.05 and c = 0.3. The continuous and 
dashed curves have the same meaning as in Fig. 1. 

Fig. 5. Suppression of the initially imposed large perturbation which is the 
fastest growing instability mode in the infinite system by the traveling pulse in 
the case a — 0.15, F — 0.2,7 = 0.05, and c = 0, in the spatial domain of the 
length L = 128 with periodic boundary conditions. The perturbation is taken 
as Upert = ^^pert = cos(fcrnaxa;) with uq — I and fcmax — 1-3. The panels (a) 
and (b) show the initial configurations of the fields, and (c) and (d) are their 
shapes produced by the direct simulations by the moment t = 420. 

Fig. 6. The instability of the pulse at the same values of parameters and for 
the same time interval as in Fig. 5, but in the spatial domain of the length four 
times as large, L — 512, without any specially imposed initial perturbation. 

Fig. 7. A stable bound state of two pulses found in the case a = 0.1, 
7 = 0.05, c = 0, r = 0.15. The panels (a) and (b) show established shapes of 
the u and v fields. 

Fig. 8. A stable bound state of three pulses in the same case as in Fig. 7. 
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